function [YinvA0,Atau,loglik,diagE,intercept,interceptE,Hlambda,E] = cicaarread(filename)
% CICAARREAD
%
%                              Y = cicaarread(filename) 
%    [invA0,Atau,loglik,diagE,intercept,interceptE,Hlambda,E] = cicaarread(filename) 
%
% inputs:
% 
%   filename : name of the file to read
%
% outputs:
% 
%          Y : matrix from file
%      invA0 : Inverse of zero-lag mixing matrix
%       Atau : Lagged mixing matrices
%     loglik : Log likelihood
%      diagE : Hessian
%  intercept : Intercept parameters
% interceptE : Hessian
%    Hlambda : Source filters
%          E : inverse hessian at A*, in packed triangular form
%
% Copyright: Mads Dyrholm, 2007
ugh = {'l','b'};
for endian = 1:length(ugh)
  [fid,msg] = fopen(filename,'r',ugh{endian});
  tmp = fread(fid,1,'int');
  if tmp==4*10
    fprintf('Endian detected.\n');
    break;
  end
  fclose(fid);
end
hdr = fread(fid,10,'int');
fread(fid,1,'int');
N = fread(fid,1,'int');
X = fread(fid,N,'real*8');
fread(fid,1,'int');
fclose(fid);
if nargout==1
  YinvA0 = reshape(X,[hdr(1) hdr(2)]);
else
  D = hdr(3);
  N = hdr(4);
  L = hdr(5);
  KC= hdr(6);
  interceptflag = hdr(8) - 1
  if interceptflag == 1
    intercept = X(D*D+D*KC*L+(1:D));
  else
    intercept = zeros(D,1);
  end
  M = hdr(9);
  version = hdr(10);
  if version<3, M=0; end
  
  dum = D*D+D*KC*L+D*M+interceptflag*D;
  YinvA0 = reshape(X(1:D*D),[D,D]);
  Atau = reshape(X(D*D+1:D*D+D*KC*L),[D,KC,L]);
  Hlambda = reshape(X(D*D+D*KC*L+1:D*D+D*KC*L+D*M),[D,M]);
  
  loglik = -X(dum+1);
  E = X(dum+2:end);
  diagE = tridiag(E,D*D+D*KC*L+D*M+D*interceptflag);
  if interceptflag==1
    interceptE = diagE(D*D+D*KC*L+D*M+(1:D));
  else
    interceptE = inf * ones(D,1);
  end
  diagE = diagE(1:D*D+D*KC*L+D*M);
end

%%%%%%%%%%%%%%5
function dig = tridiag(L,N)
inc = N;
idx = 1;
dig = zeros(1,N);
for dd = 1:N
  dig(dd) = L(idx);
  idx = idx + inc;
  inc = inc - 1;
end
